1D Haar Wavelets

This numerical tour explores 1D multiresolution analysis with Haar transform.

Installing toolboxes and setting up the path.

You need to download the following files: signal toolbox and general toolbox.

You need to unzip these toolboxes in your working directory, so that you have toolbox_signal and toolbox_general in your directory.

For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.

Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.

Execute this line only if you are using Matlab.

getd = @(p)path(p,path); % scilab users must *not* execute this

Then you can add the toolboxes to the path.

getd('toolbox_signal/');
getd('toolbox_general/');

Forward Haar Transform

The Haar transform is the simplest orthogonal wavelet transform. It is computed by iterating difference and averaging between odd and even samples of the signal.

First we load a 1D signal.

name = 'piece-regular';
n = 512;
f = rescale( load_signal(name, n) );

Initialize the transformed coefficients as the signal itself and set the initial scale as the maximum one. fw will be iteratively transformated and will contain the Haar coefficients.

fw = f;
j = log2(n)-1;

Select the sub-part of the image to transform.

A = fw(1:2^(j+1));

Compute average and differences to get coarse and details. They are weighted by 1/sqrt(2) to maintain orthogonality.

Coarse = ( A(1:2:length(A)) + A(2:2:length(A)) )/sqrt(2);
Detail = ( A(1:2:length(A)) - A(2:2:length(A)) )/sqrt(2);

Concatenate them to get the result.

A = [Coarse; Detail];

Display the result of the first step of the transform.

clf;
subplot(2,1,1);
plot(f); axis('tight'); title('Signal');
subplot(2,1,2);
plot(A); axis('tight'); title('Transformed');

Display the signal and its coarse coefficients.

s = 400; t = 40;
clf;
subplot(2,1,1);
plot(f,'.-'); axis([s-t s+t 0 1]); title('Signal (zoom)');
subplot(2,1,2);
plot(Coarse,'.-'); axis([(s-t)/2 (s+t)/2 min(A) max(A)]); title('Averages (zoom)');

Exercice 1: (the solution is exo1.m) Implement a full wavelet Haar transform that extract iteratively wavelet coefficients, by repeating these steps. Take care of choosing the correct number of steps.

exo1;

Check that the transform is orthogonal, which means that the energy of the coefficient is the same as the energy of the signal.

disp(strcat(['Energy of the signal       = ' num2str(norm(f).^2,3)]));
disp(strcat(['Energy of the coefficients = ' num2str(norm(fw).^2,3)]));

Energy of the signal       = 88.6
Energy of the coefficients = 88.6

We display the whole set of coefficients fw, with red vertical separator between the scales. Can you recognize where are the low frequencies and the high frequencies ? You can use the function plot_wavelet to help you.

clf;
plot_wavelet(fw);
axis([1 n -2 2]);

Backward Haar Transform

The backward transform computes a signal f1 from coefficients fw. If fw are the coefficients of a signal f, then one recovers f=f1. It is comuted with itertive up-sampling and filtering.

Initialize the image to recover f1 as the transformed coefficient, and select the smallezt possible scale.

f1 = fw;
j = 0;

Retrieve coarse and detail coefficients in the vertical direction.

Coarse = f1(1:2^j);
Detail = f1(2^j+1:2^(j+1));

Undo the transform by computing sums and differences.

f1(1:2:2^(j+1)) = ( Coarse + Detail )/sqrt(2);
f1(2:2:2^(j+1)) = ( Coarse - Detail )/sqrt(2);

Exercice 2: (the solution is exo2.m) Write the inverse wavelet transform that computes f1 from the coefficients fw.

exo2;

Check that we have correctly recovered the signal.

disp(strcat((['Error |f-f1|/|f| = ' num2str(norm(f-f1)/norm(f))])));

Error |f-f1|/|f| = 1.0168e-15

Haar Approximation

Non-linear approximation is obtained by thresholding low amplitude wavelet coefficients.

Set the threshold value.

T = .5;

Coefficients fw(i) smaller in magnitude than T are set to zero.

fwT = fw .* (abs(fw)>T);

Display the coefficients before and after thresholding.

clf;
subplot(2,1,1);
plot_wavelet(fw); axis('tight'); title('Original coefficients');
subplot(2,1,2);
plot_wavelet(fwT); axis('tight'); title('Thresholded coefficients');

Exercice 3: (the solution is exo3.m) Find the threshold T so that the number of remaining coefficients in fwT are a fixzd number m. Use this threshold to compute fwT and then display the corresponding approximation f1 of f. Try for an increasing number m of coeffiients.

exo3;

The Shape of a Wavelet

A wavelet coefficient fw[k] corresponds to an inner product f,psi_{j,p} where k depends on j and p.

The wavelet f = psi_{j,p} is computed by applying the inverse wavelet transform to fw where fw[k]=1 and fw[l]=0 for l \neq k.

Exercice 4: (the solution is exo4.m) Compute wavelets at several positions and scales.

exo4;